/*
 * Filename: coldplasma_main.cpp
 * Description: Simulation of electromagnetic wave propagation in a cylindrical waveguide using FDTD method.
 * 
 * 本程序实现了圆柱波导中电磁波传播的FDTD（有限差分时域）仿真。
 * 主要特点：
 * 1. 使用柱坐标系(r,θ,z)描述波导中的电磁场
 * 2. 假设场量关于θ对称，简化为二维问题
 * 3. 实现了完整的Maxwell方程组的FDTD求解
 * 4. 包含波源激励、边界条件处理和场量输出
 * 
 * 编译命令：
 * g++ -I/usr/include/eigen3 -o coldplasma_main coldplasma_main.cpp -lgsl -lgslcblas
 * 
 * 依赖库：
 * - Eigen：用于矩阵运算
 * - GSL：用于计算贝塞尔函数
 */

 #include <iostream>
 #include <Eigen/Dense>
 #include <cmath>
 #include <fstream>
 #include <gsl/gsl_sf_bessel.h> // 需安装GSL库
 
 using namespace Eigen;
 
 // 物理常数
const double c = 299792458.0;        // 光速(m/s)
const double epsilon0 = 8.854e-12;   // 真空介电常数(F/m)
const double mu0 = M_PI * 4e-7;      // 真空磁导率(H/m)
 
 struct WaveguideParams {
    // 几何参数
    double radius = 0.5;       // 波导半径(m)
    double length = 1.0;       // 轴向长度(m)
    double frequency = 2.45e9; // 频率(Hz)
    
    // 离散参数
    int Nr = 100;              // 径向网格数
    int Nz = 200;              
    
    // 计算参数
    double dr;                 // 无需默认值，由构造函数计算
    double dz;                 // 无需默认值，由构造函数计算
    double dt;                 // 由公式计算
    int num_periods = 20;      // 默认周期数
    int total_steps;           // 动态计算
    
    // 模式参数
    double kr;                 // 由贝塞尔函数根计算
    double A0 = 1.0;           // 激励幅值
    
    // 边界条件控制参数
    bool use_perfect_conductor = false; // 默认使用Mur吸收边界
    
    double kz; // 相速度参数

    // 新增等离子体参数
    double plasma_position = 0.25; // 等离子体位置 (z=L/4)
    double n_e = 0;             // 电子密度 (m^-3)

    // 精简后的初始化列表（只保留需要计算的成员）
    WaveguideParams() : 
        dr(radius/(Nr-1)),           // 先计算dr和dz
        dz(length/(Nz-1)),
        dt(1/(2*c*sqrt(1/(dr*dr)+1/(dz*dz)))), // 依赖dr/dz
        total_steps(static_cast<int>(num_periods/(frequency * dt))), // 依赖num_periods/frequency/dt
        kr(gsl_sf_bessel_zero_J1(1)/radius),    // 依赖radius
        kz(sqrt( ((2*M_PI*frequency/c)*(2*M_PI*frequency/c)) - kr*kr )) // 依赖frequency和kr
    {
        // 截止频率验证（未修改）
        double f_c = (c/(2*M_PI)) * kr;
        if(frequency < f_c) {
            std::cerr << "错误：频率低于截止频率f_c=" << f_c << "Hz，模式无法传播" << std::endl;
            exit(EXIT_FAILURE);
        }
    }

    int save_interval = 10; // 保存间隔参数
 };
 
 /**
 * @brief 场数据存储类
 * 
 * 存储波导中的电磁场分量：
 * - Etheta: θ方向电场分量
 * - Br: r方向磁感应强度分量
 * - Bz: z方向磁感应强度分量
 * 
 * 使用Eigen::MatrixXd存储二维场数据，
 * 矩阵大小为Nr×Nz，对应径向和轴向的空间离散
 */
class FieldData {
public:
    MatrixXd Etheta, Br, Bz;
    MatrixXd Jtheta; // 新增电流分量

    FieldData(int nr, int nz) : 
        Etheta(MatrixXd::Zero(nr, nz)),
        Br(MatrixXd::Zero(nr, nz)),
        Bz(MatrixXd::Zero(nr, nz)),
        Jtheta(MatrixXd::Zero(nr, nz)) {} // 初始化电流分量
};
 
 /**
 * @brief FDTD求解器类
 * 
 * 实现了FDTD算法的核心功能：
 * 1. 磁场和电场的时间步进更新
 * 2. 波源的施加
 * 3. 边界条件的处理
 * 4. 场数据的保存
 */
class FDTDEngine {
private:
    WaveguideParams params;
    
    /**
     * @brief 计算一阶贝塞尔函数J1
     * 
     * 用于计算波导中TE11模式的横向场分布
     * @param r 径向位置
     * @return J1(kr*r)的值
     */
    inline double J1(double r) { 
        return gsl_sf_bessel_J1(params.kr * r); 
    }
 
 public:
     FDTDEngine(const WaveguideParams& p) : params(p) {}
 
     /**
     * @brief 更新磁场分量
     * 
     * 基于法拉第电磁感应定律更新磁场：
     * ∂B/∂t = -∇×E
     * 
     * 在柱坐标系中展开为：
     * ∂Br/∂t = ∂Eθ/∂z
     * ∂Bz/∂t = -(1/r)∂(rEθ)/∂r
     * 
     * @param fields 场数据对象的引用
     */
    void updateMagnetic(FieldData& fields) {
        // 更新Br分量：使用中心差分格式
        for(int i=0; i<params.Nr; ++i)
            for(int j=0; j<params.Nz-1; ++j)
                fields.Br(i,j) += params.dt*(fields.Etheta(i,j+1)-fields.Etheta(i,j))/params.dz;
        
        // 更新Bz分量：包含柱坐标系的曲率项(1/r)
        for(int i=1; i<params.Nr; ++i) {
            double r = i*params.dr;
            for(int j=0; j<params.Nz; ++j)
                fields.Bz(i,j) -= params.dt*(r*fields.Etheta(i,j)-(r-params.dr)*fields.Etheta(i-1,j))/(r*params.dr);
        }
    }
 
     /**
     * @brief 更新电场分量（包含电流项）
     */
    void updateElectric(FieldData& fields, int step) {
        for(int i=1; i<params.Nr-1; ++i) {
            for(int j=1; j<params.Nz-1; ++j) {
                // 判断是否在等离子体区域内
                bool in_plasma = (j * params.dz >= params.plasma_position);

                // 更新电场分量，包含电流项
                fields.Etheta(i,j) += (params.dt/(mu0*epsilon0)) * 
                    ((fields.Br(i,j)-fields.Br(i,j-1))/params.dz - 
                    (fields.Bz(i+1,j)-fields.Bz(i,j))/params.dr);
                if (in_plasma) {
                    fields.Etheta(i,j) -= (params.dt/epsilon0) * fields.Jtheta(i,j);
                }
            }
        }

        // 在z=0处施加TE11模式波源
        double t = step * params.dt;
        for(int i=0; i<params.Nr; ++i)
            fields.Etheta(i,0) = params.A0 * J1(i*params.dr) * sin(2*M_PI*params.frequency*t);
    }

    /**
     * @brief 更新电流分量
     */
    void updateCurrent(FieldData& fields) {
        const double e = 1.602176634e-19; // 基本电荷
        const double m_e = 9.10938356e-31; // 电子质量

        for(int i=1; i<params.Nr-1; ++i) {
            for(int j=1; j<params.Nz-1; ++j) {
                // 判断是否在等离子体区域内
                bool in_plasma = (j * params.dz >= params.plasma_position);

                if (in_plasma) {
                    fields.Jtheta(i,j) += params.dt * ((e*e*params.n_e)/m_e) * fields.Etheta(i,j);
                }
            }
        }
    }

     /**
     * @brief 应用边界条件
     * 
     * 处理两类边界：
     * 1. 波导侧壁(r=a)：完美导体边界，切向电场为零
     * 2. 波导出口(z=L)：一阶Mur吸收边界条件，防止电磁波反射
     * 
     * @param fields 场数据对象的引用
     */
    void applyBoundary(FieldData& fields) {
        // 侧壁边界（r=a）：完美导体边界条件 Eθ=0
        for(int j=0; j<params.Nz; ++j)
            fields.Etheta(params.Nr-1,j) = 0.0;
        
        // 出口边界（z=L）：根据参数选择处理方式
        for(int i=0; i<params.Nr; ++i) {
            if (params.use_perfect_conductor) {
                // 理想导体边界条件：切向电场为零
                fields.Etheta(i,params.Nz-1) = 0.0;
            } else {
                // 使用相速度v_pz替代光速c
                double v_pz = (2*M_PI * params.frequency) / params.kz;
                // 计算α参数
                double alpha = (v_pz * params.dt) / params.dz;
                // 应用Mur边界条件公式
                fields.Etheta(i,params.Nz-1) = 
                    (fields.Etheta(i,params.Nz-1) + alpha * fields.Etheta(i,params.Nz-2)) 
                    / (1 + alpha);
            }
        }
    }
 
     /**
     * @brief 运行FDTD模拟（新增电流更新）
     */
    void runSimulation(FieldData& fields) {
        for(int n=0; n<params.total_steps; ++n) {
            updateMagnetic(fields);     // 更新磁场
            updateElectric(fields, n);  // 更新电场
            updateCurrent(fields);      // 更新电流
            applyBoundary(fields);      // 应用边界条件
            
            if(n % params.save_interval == 0)
                saveSnapshot(fields, n);
        }
    }
 
    void saveSnapshot(const FieldData& fields, int step) {
        // 确保data目录存在
        if (system("mkdir -p data") != 0) {
            std::cerr << "警告: 创建data目录失败" << std::endl;
        }
        
        // 保存电场数据
        char field_filename[100];
        sprintf(field_filename, "data/field_%05d.csv", step);
        std::ofstream field_out(field_filename);
        if (!field_out.is_open()) {
            std::cerr << "Failed to open file: " << field_filename << std::endl;
            return;
        }
        for (int i = 0; i < params.Nr; ++i) {
            for (int j = 0; j < params.Nz; ++j) {
                field_out << fields.Etheta(i, j);
                if (j < params.Nz - 1) {
                    field_out << ",";
                }
            }
            field_out << "\n";
        }
        field_out.close();

        // 保存电流数据
        char current_filename[100];
        sprintf(current_filename, "data/current_%05d.csv", step);
        std::ofstream current_out(current_filename);
        if (!current_out.is_open()) {
            std::cerr << "Failed to open file: " << current_filename << std::endl;
            return;
        }
        for (int i = 0; i < params.Nr; ++i) {
            for (int j = 0; j < params.Nz; ++j) {
                current_out << fields.Jtheta(i, j);
                if (j < params.Nz - 1) {
                    current_out << ",";
                }
            }
            current_out << "\n";
        }
        current_out.close();
    }
 }
 
 ;
 
 /**
 * @brief 主函数
 * 
 * 程序执行流程：
 * 1. 创建波导参数对象
 * 2. 初始化场数据（全零初值）
 * 3. 创建FDTD求解器
 * 4. 运行时域模拟
 * 
 * 运行结束后，在data/目录下生成场数据文件，
 * 可使用plot_waveguide.py脚本进行可视化
 */
int main() {
    WaveguideParams params;                    // 创建参数对象

    // 新增：计算并输出kr和kz
    double kr = params.kr; // 已由结构体初始化计算
    double omega = 2 * M_PI * params.frequency; // 角频率
    double k_z = sqrt( (omega/c)*(omega/c) - kr*kr ); // 根据用户公式计算

    std::cout << "Parameters:" << std::endl;
    std::cout << "  kr = " << kr << " 1/m" << std::endl;
    std::cout << "  kz = " << k_z << " 1/m" << std::endl;
    std::cout << "  lambda_z = " << (2 * M_PI / k_z) << " m" << std::endl; // 输出z方向波长

    // 验证传播条件（可选）
    double f_c = (c/(2*M_PI)) * kr; // 截止频率
    if(params.frequency < f_c) {
        std::cerr << "警告：频率低于截止频率f_c=" << f_c << "Hz，可能导致无解" << std::endl;
    }

    // 新增参数设置示例（可选）
    params.use_perfect_conductor = false; // 是否开启理想导体边界
    
    // 输出使用的边界条件
    if (params.use_perfect_conductor) {
        std::cout << "使用边界条件: 理想导体" << std::endl;
    } else {
        std::cout << "使用边界条件: Mur吸收边界" << std::endl;
    }
    
    FieldData fields(params.Nr, params.Nz);    // 初始化场数据
    FDTDEngine simulator(params);              // 创建FDTD求解器
    
    simulator.runSimulation(fields);           // 运行模拟
    return 0;
}
 
 
 
 